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FORCE EVALUATION IN THE LATTICE BOLTZMANN METHOD INVOLVING 

CURVED GEOMETRY 

RENWEI MET, DAZHI YU*, WEI SHYY* . AND LI-SHI LUO§ 


Abstract. The present work investigates two approaches for force evaluation in the lattice Boltzmann 
equation: the momentum-exchange method and the stress-integration method on the surface of a body. 
The boundary condition for the particle distribution functions on curved geometries is handled with second 
order accuracy based on our recent works. The stress-integration method is computationally laborious for 
two-dimensional flows and in general difficult to implement for three-dimensional flows, while the momentum- 
exchange method is reliable, accurate, and easy to implement for both two-dimensional and three-dimensional 
flows. Several test cases are selected to evaluate the present methods, including: (i) two-dimensional pressure- 
driven channel flow; (ii) two-dimensional uniform flow past a column of cylinders; (iii) two-dimensional flow 
past a cylinder asymmetrically placed in a channel (with vortex shedding); (iv) three-dimensional pressure- 
driven flow in a circular pipe; and (v) three-dimensional flow past a sphere. The drag evaluated by using 
the momentum-exchange method agrees well with the exact or other published results. 

Key words, lattice Boltzmann method, force evaluation on fluid-solid interface, momentum-exchange 
method, stress-integration method, boundary condition for curved geometries, accuracy, 3-D flows 

Subject classification. Fluid Mechanics 
1. Introduction. 

1.1. Background of the lattice Boltzmann equation method. The method of lattice Boltzmann 
equation (LBE) solves the microscopic kinetic equation for particle distribution function /(«,£,<), where £ 
is the particle velocity, in phase space (as, £) and time /, from which the macroscopic quantities (flow mass 
density p and velocity «) are obtained through moment integration of Because the solution pro- 

cedure is explicit, easy to implement and parallelize, the LBE method has increasingly become an attractive 
alternative computational method for solving fluid dynamics problems in various systems [1, 2, 3, 4]. The 
most widely used lattice Boltzmann equation [1, 2, 3, 4] is a discretized version of the model Boltzmann 
equation with a single relaxation time approximation due to Bhatnagar, Gross, and Krook (BGK model) 

[5], 

^/ + £V/ = i[/-/W], (1.1) 

where /'°' is the Maxwell-Boltzmann equilibrium distribution function and A is the relaxation time. The 
mass density p and momentum density pu are the first. (D + 1) hydrodynamic moments of the distribution 
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function / and / !0 > in D dimensions. It can be shown that the particle velocity space £ can be discretized 
and reduced to a very small set of discrete velocities {£<*|a = 1,2, and the hydrodynamic moments 

of / and /'°) as well as their fluxes can be preserved exactly, because the moment integral can be replaced 
by quadrature exactly up to a certain order in £ [6, 7, 8, 9]. With velocity space £ properly discretized, 
Eq. (1.1) reduces to a discrete velocity model of the Boltzmann equation 

$/a+€a-V/ a = i [/„-/<?>], (1.2) 

In the above equation, f a (x,t) = /(#,£«, t) and fa\x,t) = f ( 0 \x,£ a ,t) are the distribution function and 
the equilibrium distribution function of the ath discrete velocity £ a , respectively. Equation (1.2) is then 
discretized in space x and time t into 

fjxi + e a S U t + St) - f a (Xi , t) - f)] . (1.3) 

T 

where r = X/St is the dimensionless relaxation time and e a is a discrete velocity vector. The coherent 
discretization of space and time is done in such a way that Sx = e a St is always the displacement vector from 
a lattice site to one of its neighboring sites. The equilibrium distribution function fa q \xi,t) in the lattice 
Boltzmann equation (1.3) is obtained by expanding the Maxwell-Boltzmann distribution function in Taylor 
series of u up to second order [6, 7], and can be expressed in general as 


f(eq) 
J a 


= W a p 


1 + ^(e a -u) + ~(e a ■ u f 


3 

2c 2 



(1.4) 


where c = S x /St ; S x is the lattice constant of the underlying lattice space; and coefficient w a depends on the 
discrete velocity set {e a j in D spatial dimensions. In what follows, we shall use the lattice units of S x = 1 
and St = 1. The Appendix provides the det ails of coefficient w a and the discrete velocity set {e tt } for the two- 
dimensional nine-velocity model (D2Q9) and the three-dimensional nineteen-velocity model (D3Q19) [10]. 
Figure 1 shows the discrete velocity sets of the two models. It should be pointed out that there exist other 
discrete velocity sets {e a \ that have the sufficient symmetry for hydrodynamics [6, 7]. A comparative study 
of three three-dimensional LBE models including the fifteen-velocity model (D3Q15), the nineteen-velocity 
model (D3Q19), and the twenty-seven-velocity model (D3Q27), in terms of accuracy and computational 
efficiency has been conducted by Mei et al. [ 11 ]. It was found that the nineteen-velocity model (D3Q19) 
offers a better combination of computational stability and accuracy. The D2Q9 and D3Q19 models will be 
used in this study for force evaluation in two-dimensional arid three-dimensional flows, respectively. Equation 
(1.3) is conveniently solved in two steps 

collision: f a (xi, t ) = /„(aq, t) - - [/«(»*, t) - /i eq) (aq, <)] , (1.5a) 

streaming: f a (xi + e a S t , t + S t ) = f a (xi, t) , (1 .5b) 


which is known as the LBGK scheme [1, 2]. The collision step is completely local and the streaming step 
is uniform and requires little computational effort, which makes Eq. (1.5) ideal for parallel implementation. 
The simplicity and compact nature of the LBGK scheme, however, necessitate the use of the square lattices 
of constant spacing (<5 X = S y ). and consequently lead to the unity of the local Courant-Friedrichs-Lewy 
number, because St = 5 X = 1. 

1.2. Boundary condition for a curved geometry in the LBE method. Consider a part of an 
arbitrary curved wall geometry, as shown in Fig. 2, where the filled small circles on the boundary, x w , denote 
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Fig. 1. Discrete velocity set {e a }. (left) Two-dimensional nine-velocity (D2Q9) model, (right) Three-dimensional 
nineteen-velocity (D3Q19) model. 


the intersections of the boundary with various lattice-to-lattice links. The fraction of an intersected link in 
the fluid region, A, is defined by 


^ \\ x ,f fly 1 1 

II*/ - *tll ' 


( 1 . 6 ) 


Obviously the horizontal or vertical distance between an, and x w is A5 C on the square lattice, and 0 < A < 1. 
In Eq. (1.5b), the value of f a (xi,t) needs to be constructed according to the location of the boundary and 
the boundary conditions, if the grid point aq = an, lies beyond the boundary. In the past, the bounce- 
back boundary condition has been use to deal with a solid boundary in order to approximate the no-slip 
boundary condition at the solid boundary [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. However, it is 
well understood that this bounce-back boundary condition satisfies (die no-slip boundary condition with a 
second-order accuracy (for the Couette and Poiseuille flows) at the location one half lattice spacing (A = 1 /2) 
outside of a boundary node where the bounce-back collision takes place; and this is only true with simple 
boundaries of straight lines parallel to the lattice grid [19, 20, 21]. For a curved geometry, simply placing the 
boundary halfway between two nodes wall alter the geometry on the grid level and degrade the accuracy of 
the flow field and the force on the body at finite and higher Reynolds number. To circumvent this difficulty, 
Mei and Shyy solved Eq. (1.2) in curvilinear coordinates using a finite difference method to compute f a 
[25]. He and Doolen used body-fitted curvilinear coordinates with interpolation throughout the entire mesh, 
except at the boundaries where the bounce-back boundary condition is used [26]. In the recent works of 
Filippova and Han el [27] and Mei et al. [28, 11], a second-order accurate boundary condition for curved 
geometry was developed in conjunction with the use of Cartesian grids in order to retain the advantages 
of the LBE method. An interpolation scheme is employed only at the boundaries to obtain f a (xi,t). The 
detailed assessment on the impact of the boundary condition on the accuracy of the flow' field has been given 
in Ref. [28] for some two-dimensional flows and in Ref. [11] for some three-dimensional flow's. 

Because the bounce-back type boundary conditions play an important role iri lattice Boltzmann simula- 
tions, it is important for us to understand how the boundary conditions work. First of all, one must realize 
that it is impossible for any kinetic numerical scheme to impose a given velocity (the Dirichlet boundary- 
condition) on a given grid node, because the Knudsen layer type of phenomena [29, 30, 31] would be mani- 
fested in kinetic schemes [32, 19, 20, 21]. For example, in the Poiseuille and the Couette flows, the location 
where hydrodynamic boundary conditions are satisfied are one-half grid spacing away from the boundary- 
grids where the bounce-back boundary conditions are imposed [19, 20, 21]. For flows around an arbitrary- 
shaped body analytical solutions do not exist. Nevertheless, substantial evidence show's that the bounce- 
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.Fig. 2. layout of the regularly spaced lattices and curved wall boundary. The circles (o), discs (•), shaded discs («), and 
diamonds (O) denote fluid nodes, boundary locations (*«.), solid nodes which are also boundary nodes (*;,) inside solid, and 
solid nodes, respectively. 


back boundary conditions combined with interpolations, and including the one-half grid spacing correction 
at boundaries, are in fact second-order accurate and capable of handling curved boundaries [26, 23, 24, 33]. 
This point is also demonstrated in the present work. 


1.3. Force evaluation and related works. In spite of numerous improvement for the LBE method 
during the last several years, one important issue that has not, been systematically studied is the accurate 
determination of the fluid dynamic force involving curved boundaries. Needless to say, accurate evaluation 
of the force is crucial to the study of fluid dynamics, especially in fluid-structure interaction. Several force 
evaluation schemes, including momentum exchange [14, 16] and integration of surface stress [26, 34], have 
been used to evaluate the fluid dynamic force on a curved body in the context of the LBE method. 

He and Doolen [26] evaluated the force by integrating the total stress on the surface of the cylinder and 
the components of the stress tensor were obtained by taking respective velocity gradients. Even though a 
body-fitted grid was used, an extrapolation was needed to obtain the stress in order to correct the half- 
grid-cell spacing effect due to the bounce-back boundary condition. Filippova and Hanel [27] developed 
a second-order accurate boundary condition for curved boundaries. However, the fluid dynamics force on 
a circular cylinder asymmetrically placed in a. two-dimensional channel was obtained by integrating the 
pressure and deviatoric stresses on the surface of the cylinder by extrapolating from the nearby Cartesian 
grids to the solid boundary [27, 34]. To gain insight into the method of surface stress integration, it is 
instructive to examine the variation of the pressure on the surface of a circular cylinder at finite Reynolds 
number obtained by using the LBE method for flow over a column of cylinders (see Ref. [28], and Sec. 3.2). 
Figure 3 shows the pressure coefficient, 


Cp 


P ~ /'.. 
§ pU' 2 
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Fig. 3. Flow past a column of 2D circular cylinders. Distribution of the pressure coefficient Cp on the surface of a 2D 
circular cylinder of radius r — 6.6, and center-to- center distance H/r — 20. The stagnation point, is located at 0 — 180° . The 
LBE result denoted by symbols x is obtained with r — 0.6 and Re — 40. The solid line is the result obtained by using a 3D 
multi-block, body-fitted grid, and pressure-based Navier-Stokes solver with a much finer resolution. 


on the surface obtained by using second-order extrapolation, where p 0 0 is the far upstream pressure. Only 
those boundary points, x w , intersected by the horizontal or vertical velocities, i.e., e\, e 3 , eg, and 67 , are 
considered in the result given by Fig. 3. If the boundary points intersected by the links in the diagonal 
velocities, i.e., e. 2 , e*, eg, and eg, are also considered, the variation of Cp would be more noisy. The 
components of the deviatoric stress tensor show a similar noisy pattern. It is not clear how the noise in 
the pressure and stresses affect the accuracy of the fluid dynamic force in the stress-integration method. 
While the programming in the extrapolation and integration is manageable in two-dimensional (3D) cases, 
it is rather laborious in three-dimensional cases. In Fig. 3, the LBE result of Cp(0) (indicated by symbol 
x) is compared with that obtained by using a 3D multiblock, body-fitted coordinates, and pressure-based 
Navier-Stokes solver [35, 36, 37] with a much finer resolution: 201 points around the cylinder and the smallest 
grid size along the radial direction dr — 0.026 (relative to r = 1). Not surprisingly, the result obtained by 
using the Navier-Stokes solver with body-fitted grid and high resolution is smoother than the LBE result 
with a Cartesian grid of coarser resolution. Nevertheless, the LBE solution still essentially agrees with the 
N avier- Stokes solution. 

Instead of the stress-integration method, Ladd used the momentum-exchange method to compute the 
fluid force on a sphere in suspension flow [14], In the flow simulation using the bounce-back boundary 
condition, the body is effectively replaced by a series of stairs. Each segment on the surface has an area of 
unity for a cubic lattice. The force on each link [halfway between two lattices at x / and xy = (xf + e a St) 
in which xy resides in the solid region] results from the momentum exchange (per unit time) between two 
opposing directions of the neighboring lattices 

7 -[e a fa(xf) - eafa(xf + e a St)] 

<Jt 

in which e a = — e a . Whereas the momentum-exchange method is very easy to implement computationally, 
its applicability and accuracy for a curved boundary have not been systematically studied. To recapitulate, 




there are two major problems associated with the method of surface stress integration. First, the components 
of stress tensor are often noisy on a curved surface due to limited resolution near the body and the use of 
Cartesian grids. The accuracy of such a method has not been addressed in the literature. Second, the 
implementation of the extrapolation for Cartesian components of the stress tensor to the boundary surface 
and the integration of the stresses on the surface of a three-dimensional geometry are very laborious in 
comparison with the intrinsic simplicity of the lattice Boltzmann simulations for flow' field. The problems 
associated with the method of momentum exchange are as follows, (a) The scheme was proposed for the case 
with A = 1/2 at every boundary intersection x w . Whether this scheme can be applied to the cases where 
A yf 1/2, when, for example, the boundary is not straight, needs to be investigated, (b) As in the case of 
stress-integration method, the resolution near a solid body is often limited and the near wall flow' variables 
can be noisy. If one uses the momentum-exchange method to compute the total force, it is not clear what 
the adequate resolution is to obtain reliable fluid dynamic force on a bluff body at a given (moderate) value 
of Reynolds number, say. Re « O(10 2 ). 

1.4. Scope of the present work. In what follows, two methods for the force evaluation, i.e. , the 
stress-integration and the momentum-exchange methods, will be described in detail. The shear and normal 
stresses on the wall in a pressure driven channel flow will be first examined to assess the suitability of the 
momentum-exchange method when A / 1/2 and analyze the errors incurred. The results on the drag force 
for flow' over a column of circular cylinders using these two methods will be subsequently assessed for the 
consistency. The drag coefficient at Re = 100 are compared with the result of Fornberg [38] obtained by 
using a second-order accurate finite difference scheme with sufficient grid resolution. For flow over a cylinder 
asymmetrically placed in a channel at Re = 100, the unsteady drag and lift coefficients are computed and 
compared with the results in the literature. The momentum-exchange method is further evaluated for three- 
dimensional fully developed pipe flow' and for a uniform flow over a two-dimensional array of spheres at finite 
Reynolds number. We found that the simple momentum-exchange method for force evaluation gives fairly 
reliable results for the two-dimensional and three-dimensional flows. 

2. Methods for Force Evaluation in LBE Method. 


2 .1. Second-order accurate no-slip boundary condition for curved geometry. The analysis of 
boundary conditions for a curved boundary in the lattice Boltzmann equation is accomplished by applying 
the Chapman-Enskog expansion for the distribution function at the boundary. The following approximation 
for the post-collision distribution function on the right-hand side of Eq. (1.5b) can lead to a second-order 
accurate no-slip boundary condition [11, 27, 28] 

_ 3 

fa(xb, t ) = (1 - x)f a (x f , t) + xfa( x b , t) + 2w a p~ea ■ U w , (2.1) 

where 


fa(XbJ-) = W a p(Xf,t) 


3 9 3 

1 + ^(e a • u bf ) + ^(e a • Uff - 


= fa q) (xf, t) + w a p(x f , t) —e a • (u b f - Uf) , 


( 2 . 2 ) 


and 


u b f = Uff = Uf(x f + e & S t ,t ) , 


(2A - 1) 
(r-2) 


1 3 

Ubf = 7m~( 2A - 3 )u f + —u w , 


(2 A - 1) 

(V + 1/2) 


0<A<-, 


- < A < 1. 


(2.3a) 

(2.3b) 
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The above treatment is applicable for both the two-dimensional and three-dimensional lattice Boltzmann 
models. 


By substitution of Eq. (2.2), Eq. (2.1) becomes 


f s (Xb, t) = f a (Xf,t) ~ X /«(*/, t) - f ( a eq) (Xf, t) 


+w a p(xf, t)-^e a ■ ( u b f -Uf- 2u w ). 


(2.4) 


Thus, the above treatment of curved boundaries can be thought as a modification of the relaxation (the 
viscous effect) near the wall (with the relaxation parameter x), in additional to a forcing term accounting 
for the momentum-exchange effect due to the wall. 


2.2. Force evaluation based on stress integration. He and Doolen [26] evaluated the force by 
integrating the total stresses on the boundary of the cylinder d fl, 

F=f d.4w-{-pl + p/v[(V:u) + (V:ii) T ]} , (2.5) 

Jon 

where n is the unit out normal vector of the boundary dfl, I is the identity tensor of second rank, V : u 
denotes the second rank tensor whose components are and T is the tranpose operator. In Ref. [26], a 
body-fitted coordinate system together with grid stretching was used such that a large number of grids can 
be placed near the body to yield a reliable velocity gradient d In general, since u is not the primary 
variable in the LBE simulations and the evaluation of u using e a f a based on f a ’s suffers the loss of 
accuracy due to the cancellation of two close numbers in /«’ s the evaluation of the derivative d,uj will 
result in further degradation of the accuracy. Filippova [34] used a similar integration scheme to obtain the 
dynamic force on the body for the force on a circular cylinder [271 except that the deviatoric stresses were 
evaluated using the non-equilibrium part of the particle distribution function [see Eq. (2.7) below]. However, 
since a Cartesian grid was used, the stress vectors on the surface of the body (with arbitrary A) have to 
be computed through an extrapolation procedure based upon the information in the flow field. This leads 
to further loss of accuracy for a finite lattice size S x when the shear-layer near the wall is not sufficiently 
resolved. 

In Eq. (2.5), the pressure p can be easily evaluated using the equation of state p = <? a p. For D2Q9 and 
D3Q19 models, c% = 1/3 so that p — p/3. The deviatoric stress for two-dimensional incompressible flow 


nj = pv ( diUj + djUi) (2.6) 

can be evaluated using the non-equilibrium part of the distribution function /£ n * q ! — [/ Q — /)® q) ] 

Tij = f l /« neq) (*> *) ■ e a SijJ , (2.7) 

where e a ,i and e a j are ith and jth Cartesian component of the discrete velocity e a , respectively. For the 
flow past a circular cylinder, a separate set of surface points on the cylinder can be introduced in order 
to carry out the numerical integration given by Eq. (2.5). The values of the pressure and each of the six 
components of the symmetric deviatoric stress tensor on the surface points can be obtained using a second- 
order extrapolation scheme based on the values of p and nj at the neighboring fluid lattices. The force 
exerting on the boundary Oil is computed as 

F = / d-Aw • {—pi + pl/[(V:rt) + ( V : xt) T ] } extrapo i ate d • (2.8) 

J oil 

It is worth commenting here that for the two-dimensional flow past a cylinder, nearly half of the length of 
the entire code is taken up by the above force evaluation procedure. 



2.3. Method based on the momentum exchange. In order to employ the momentum-exchange 
method efficiently, two scalar arrays, w(i,j) and w b (i,j) are introduced. A value of 0 is assigned to w(z,j) 
for the lattice site (i,j) that are occupied by fluid; a value of 1 is assigned to w(i,j) for those lattice nodes 
inside the solid body. The array is set to zero everywhere except for those boundary nodes, an,, where 

a value of 1 is assigned. For a given nonzero velocity e a , e s denotes the velocity in the opposite direction, 
i.e. , e s = —e a (see Fig. 2). For a given boundary node x b inside the solid region with w b (i,j) = 1 and 
w(i,j) = 1, the momentum exchange with all possible neighboring fluid nodes over a time step St = 1 is 

Y2 e a [fa(xb, t ) + fa{x h + e & 6 t , t ) j [1 -- w(x b + e a S t )] ■ 

a^O 

Simply summing the contribution over all boundary nodes x b belonging to the body, the total force (acted 
by the solid body on the fluid) is obtained as 

F = YL Y2 ea [/«(*»>*) + /»(*/, + e a <M)] [t- w(x b + e a d'<)]. (2.9) 

all a^O 

In the momentum-exchange method the force F is evaluated after the collision step is carried out and 
the value of / a at the boundary given by Eq. (2.1) has been evaluated. The momentum exchange occurs 
during the subsequent streaming step when f s (x b ,t + St) and f a (xf, t + St) move to x / and x b , respectively. 
As mentioned in the introductory section, the effect of the variable A is not explicitly included, but it is 
implicitly taken into account in the determination of f & (x b , t + St). The applicability of Eq. (2.9) will be 
examined and validated. 

Clearly, the force is proportional to the number of boundary nodes x b in the above formula of F and 
the number of the boundary nodes increase linearly with the size of the body in a two-dimensional flow. 
However, since the force is normalized by pU 2 r in the formula for Cd in two-dimensions [see Eq. (3.9)], the 
drag coefficient Cd should be independent of r. 

3. Results and Discussions. For straight walls, there is no doubt that Eq. (2.5) together with the 
equation of state for pressure and Eq. (2.7) for ry gives accurate results for the force provided that the /„’ s 
are accurately computed. To demonstrate the correctness of Eq. (2.9) based on the momentum exchange for 
an arbitrary A, we first consider the pressure driven channel flow (see Fig. 4) for which exact solutions for 
the velocity and stresses are known. The second case considered is the two-dimensional flow past a column of 
circular cylinders at Reynolds number Re = 100 and H/r = 20, where H is the distance between the centers 
of two adjacent cylinders. The values of the drag computed using the two force evaluation methods are then 
compared with the result of Fomberg [38]. The dependence of the drag on the radius r in the momentum- 
exchange method is examined to assess the reliability of this method. The third case is the two-dimensional 
flow over a circular cylinder that is asymmetrically placed in a channel at Re = 100 (with vortex shedding). 
The time dependence of the drag and lift coefficients is compared with results in the literature. 

We also consider two cases of three-dimensional flow. The first case is the pressure driven flow in a 
circular pipe for which the exact solutions for both the velocity profile and the wall shear stresses are known. 
The assessment for the momentum- exchange method for three-dimensional flows will be made first in this 
case. Finally, the momentum-exchange method will be evaluated by considering the drag on a sphere due to 
a uniform flow over a sphere in a finite domain. The details for the flow field computation can be found in 
Refs, [28. 11]. 

3.1. Two-dimensional pressure-driven channel flow. In the case of the channel flow, the force 
on the top wall (y = H ) at a given location x (i = N x /2 + 1, for example) can be evaluated using the 
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Fig. 4. The channel flow configuration in the LBE simulations with an arbitrary A. 


momentum-exchange method as follows. The wall is located between j = N y and N y — 1 (cf. Fig. 4). The 
x and y components of the force on the fluid at the top wall near the ith node are 

F x = [Mi, 3) - /a(* - 1,3 ~ 1)] e 9 , x + [fs(i,j) + fi{i + 1,3 ~ 1)] e 8 , x (3.1a) 

F y = [Mi, 3) + h(i - 1 ,3 ~ 1)] e s, y + [h (i,j) + h(i + 1 ,j - 1)] e»,„ 

+[h(i,j) + Mi, 3 ~ 1)] er, y , (3.1b) 


where c a j denotes the jth Cartesian component of velocity e a . Since 6 X = 1, F x and F v are, effectively, the 
total shear and normal stresses, a xy and a yy , which include the pressure and the deviatoric stresses, on the 
fluid element at y = H. 

Based on Eq. (2.7), the deviatoric component of the fluid shear stresses at j = N y - 1 (or y — N v — 3 + A) 
and N y - 2 (or y — N y - 4 + A) can be exactly evaluated based on the nonequilibrium part of the distribution 
functions in the flow field if they are correctly given. A linear extrapolation of the deviatoric shear stresses 
to y = H = N y — 3 + 2A yields 

- T X y(j - Ny - 1) + A [T Xy {j - Ny ~ l) ~ T X y (j = Ny ~ 2)] , (3.2) 


where the superscript “(neq)” denotes the value computed from /a” eq , and the subscript w refers to the 
value at the wall. The deviatoric normal stress, T yy *w, can be similarly computed. In a fully developed 
channel flow, the normal component of the deviatoric stress T vy {y) is expected to be zero while the total 
normal stress <r yy (y) is equal to the negative of the pressure (— p). It needs to be pointed out that this 
method of evaluating r^w given by Eq. (3.2) for two-dimensional channel flow is equivalent to the method 
of the surface stress integration based on the extrapolated pressure and the deviatoric stresses on the solid 
wall except that no numerical integration on the solid surface is needed. 

After the velocity profile u x (y) is obtained from f a , the shear stress r xy on the wall can also be calculated 
using the near wall velocity profile as 


(2 + A) [0 - u x {j = N y - 1)] 
Pi/ (1 + A) A 

A 


(1 + A) 


[«x(j = Ny - 1) - U X (j = Ny - 2)] . 


In the above, a linear extrapolation is employed to evaluate the velocity derivative (du x /dy)\ y= n at the wall. 
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Table 3.1 


Comparison of fluid stresses at y — H in a two-dimensional pressure driven channel flow with dp/dx — —1.0 x 10 6 in 
the lattice units, N y — 35 and r — 0.6 as a function of A. Column 2. given by Eq. (3.4); Column 3, —F x given 

by Eq. (3.1a); Column 4 ? ~ T iy^w' given by Eq. (3.2); Column 5, —po(du x fdy)\ v —h Eq. (3.3); Column 6, —F y given by 
Eq. (3.1b); Column 7, pressure p obtained in the simulation. 


A 

exact v i no 
~ T xy,w x 

-F x x 10 5 

neq w in5 
T xy,w x 

- \ v // x 

-K 

P 

0.01 

1.601 

1.6333 

1.6010 

3.5294 

0.3333 

0.3333 

0.02 

1.602 

1.6333 

1.6020 

2.5555 

0.3333 

0.3333 

0.03 

1.603 

1.6333 

1.6030 

2.2309 

0.3333 

0.3333 

0.04 

1.604 

1.6333 

1.6040 

2.0685 

0.3333 

0.3333 

0.05 

1.605 

1.6333 

1.6050 

1.9710 

0.3333 

0.3333 

0.1 

1.610 

1.6333 

1.6100 

1.7760 

0.3333 

0.3333 

0.2 

1.620 

1.6333 

1.6200 

1.6781 

0.3333 

0.3333 

0.25 

1.625 

1.6333 

1.6250 

1.6583 

0.3333 

0.3333 

0.3 

1.630 

1.6333 

1.6300 

1.6451 

0.3333 

0.3333 

0.3333 

1.633 

1.6333 

1.6330 

1.6385 

0.3333 

0.3333 

0.35 

1.635 

1.6333 

1.6350 

1.6357 

0.3333 

0.3333 

0.4 

1.640 

1.6333 

1.6400 

1.6285 

0.3333 

0.3333 

0.5 

1.650 

1.6333 

1.6500 

1.6184 

0.3333 

0.3333 

0.6 

1.660 

1.6333 

1.6600 

1.6214 

0.3333 

0.3333 

0.7 

1.670 

1.6333 

1.6700 

1.6244 

0.3333 

0.3333 

0.8 

1.680 

1.6333 

1.6800 

1.6274 

0.3333 

0.3333 

0.9 

1.690 

1.6333 

1.6900 

1.6305 

0.3333 

0.3333 

0.95 

1.695 

1.6333 

1.6950 

1.6321 

0.3333 

0.3333 

0.99 

1.699 

1.6333 

1.6990 

1.6335 

0.3333 

0.3333 


Finally, the exact solution for the fluid shear stress on the wall (y — H) is 

= \% H > H = Ny ~ 3 + 2A (3.4) 

based on the parabolic velocity profile or simple control volume analysis. This exact result can be used to 
assess the accuracy of the aforementioned methods for the force evaluation. 

In the LBE simulations, the pressure gradient is enforced through the addition of an equivalent body 
force after the collision step [26, 11]. While the velocity field given by the LBE solution can be unique, 
the pressure field [thus the density field p(x,y)\ can only be unique up to an arbitrary constant. In view 
of Eq. (3.3), it is difficult to compare the stresses for different cases if p(i,j ) converges to different values 
in each case. To circumvent this difficulty, the density field in the channel flow simulation is normalized by 
pit — 2,j = Ny/2) at every time step. This normalization procedure results in p(x,y ) = 1 throughout the 
entire computational domain. It is also applied to the three-dimensional flow in a circular pipe. 

Table 3.1 compares the numerical values of the shear stress for a typical case (N y = 35, dpjdx — -TQ~ 6 
in the lattice units, and r = 0.6) based on: given by Eq. (3.4), F x given by Eq. (3.1a), given 

by Eq. (3.2), and pv{du x j dy)\ v= n given by Eq. (3.3). Also listed is the comparison between F y given by 
Eq. (3.1b) and —p. All computations are carried out with double precision accuracy. 

It is noted that riy^ is identical to for all values of A. A closer examination of the shear 
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stress profile using Eq. (2.7) across the channel reveals that T^wiy) is also equal to the exact shear stress 
profile r“ act (y), which is linear, despite the errors in the velocity profile u x (y) for all values of A. A linear 
extrapolation, Eq. (3.2), for a linear profile therefore gives the exact wall shear stress. Thus, the exactness 
of Txy^.w i n the LBE simulation of channel flow indicates the reliability of the LBE solution for the stress 
field r-j eq) (x,y) by using Eq. (2.7). However, as Fig. 3 indicates, the accuracy of integrating to 

obtain the fluid dynamic force in nontrivial geometries needs to be further investigated, as will be discussed 
in the following sections. 

For 0 < A < 1, the normal force F y given by Eq. (3.1b) based on the momentum- exchange method 
agrees exactly with the pressure on the wall. This is a rather special quantity since deviatoric component of 
the force is identically zero. Nevertheless, the method of the momentum exchange does give a reliable value 
for the normal stress. 

For the shear (tangential) force, it is observed from Table 3.1 that for fixed dp/dx, F x does riot change 
as A increases from 0.01 to 0.99. On the other hand, the exact result = | (dp/dx)(N y — 3 + 2A), 

increases linearly with A. Further computations were carried out over a range of N y (= 35, 67, 99, and 
131) and r (= 0.505, 0.51, 0.52, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, and 1.6). The results indicate that the 
momentum-exchange method gives the shear stress on the top wall as 


F x 


1 dp 

2 dx 



(3.5) 


That is, F x is independent of r and A. The error in F x is zero when A = 1/3. The absolute error attains 
the maximum when A = 1, which gives the relative error of 4/3 H for F x . Although the frequently used 
momentum-exchange method is a natural choice for the force evaluation in conjunction with the bounce-back 
boundary condition for A = 1/2, one must be aware that this method is not exact and the error in the force 
evaluation using the momentum- exchange method depends on A and the resolution. 

The error in F x is due to the fact that the derivatives of the velocity field are not considered in the 
boundary conditions. This can be understood by analyzing Eq. (3.1a). At the steady state, and with the 
approximation that 

fa « /£**’ + fa ] = /i eq) - ^WoP : ( |(e« • V)(e« - u ) , (3.6) 

Equation (3.1a) at the top wall becomes 

3 

F x ss 2u’2p-je2 ■ (ut,f + Uf — 2u w ) , (3.7) 

where the substitution of Eq. (2.4) for / 6 and / 8 has been made. The only term in the above equation which 
has A dependence is Uf,f. When 0 < A 1/2, F x is independent of A, and when 1/2 < A < 1, F x weakly 
depends on A because u w = 0 in this case [see Eqs. (2.3)]. In the case where F x is obtained by summing 
over a set of symmetric lattice points, cancellations in the summation may further weaken the dependence 
of F x on A. 

Table 3.1 also shows that for the shear stress based on the derivative of the velocity obtained by using 
finite-difference, the loss of accuracy is quite significant for small values of A (< 0.05) when r = 0.6. For 
other values of A (> 0.3), the accuracy is comparable with that of F x . However, as shown in Fig. 5(a), 
the accuracy of pv{du x / dy)\ v= n based on the near- wall velocity derivative deteriorates as the relaxation 
time r increases (from 0.51 to 1.6). To see the cause of the increasing error in pv(du x j dy)\ v= n , Fig. 5(b) 
shows dimensionless wall velocity, u w /u c , obtained by a three-point second-order Lagrangian extrapolation 
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Fig. 5. The LBE simulations of the channel flow . with A = 0.2, 1/3, 0.5, and 0.7. The pressure drop is d x p — — l.Ox 10"' 6 
in lattice units, (a) Ratio between the wall force, pvd v u x \ y =H , evaluated by using Eg. (3.3), and the exact value = 

—Hdxp/2, given by Eg. (3.4) as a Junction of r. (b) Normalized wall slip velocity u w /u c as a function of r. 


of the near wall velocity u x (y) as a function of r. The increasing slip velocity u w on the wall with the 
increasing relaxation time r was also observed in Ref. [15]. It is the result of increasing particle mean free 
path that causes the deviation of the kinetic solution from the hydrodynamic solution. It is clear that the 
poor performance of pv{du x / dy)\ v =n is associated with the increasing error in the near wall velocity profile 
as t increases. Since the stress tensor ry can be calculated directly from f a [see Eq. (2.7)] without the 
need for directly computing velocity derivatives, the force evaluation method based on the evaluation of the 
velocity gradient in the form of Eq. (2.6) is not recommended. 


3.2. Steady uniform flow over a column of cylinders. For a uniform flow over a. column of circular 
cylinders of radius r and center-to-center distance H (see the left part of Fig. 6 for illustration), symmetry 
conditions for f a ' s are imposed at y — ±H/2. Most of the details of flow field simulation can be found in 
Ref. [28]. The Reynolds number is defined by the diameter of the cylinder d as Re — Ud/v , where U is the 
uniform velocity in the inlet. It must be noted that for a consistent determination of the force, the upstream 
boundary must be placed far upstream. A shorter distance between the cylinder and the boundary will result 
in higher drag. In this study, it is placed at about, 20 radii to the left of the center of the cylinder. Reducing 
the distance between the boundary and the cylinder to 12.5 radii while keeping the rest of the computational 
parameters fixed would increase the drag coefficient by about 1.8% at Re = 100. The downstream boundary 
is located about 25 - 30 radii behind the cylinder to allow sufficient wake development. The simulation is 
terminated when the following criterion based on the relative 1 , 2 -norm error in the fluid region O is satisfied, 


Eo = 


yi \\u{xi,t + 1 ) — u(xi,t)\\ 2 

aj;6£2 


N 


\\u{Xi,t + 1)|| 2 

X{ GQ 


< e. 


(3.8) 


In this case, e = 10 6 was chosen for both Re = 10 and 100. 

Following Fornberg [38], the drag coefficient over a circular cylinder of radius r is defined as 


C-D = 


m 

pU 2 r 


(3.9) 


Figure 7(a) compares C-d obtained from: momentum-exchange method, surface stress integration, and finite 
difference result of Fornberg [38] using a vorticity-stream function formulation at Re = 100, H/r = 20, and 
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Fig. 6. Computational domain for the uniform flow past a sphere of radius r. The dashed lines indicate boundaries of 
computational domain, (left) Unbounded domain in xy plane, and (right) bounded domain in yz plane. 

radius r ranging from 2.8 to 13.2. For r > 8, both the momentum exchange and the stress integration 
methods give satisfactory results for Cd in comparison with the value of 1.248 given in Ref. [38]. The 
small differences in Cd could be due to the fact that in Ref. [38], the computational domain is much larger 
in the downstream direction the downstream boundary condition is imposed at 300 radii behind the 
cylinder in Ref. [38], as opposed to 25 - 30 radii here. This adds credence to the validity of Eq. (2.9) for 
evaluating the total force on a body. The values of Cd from the momentum-exchange method have a little 
less variation than that from the stress integration. Accepting an error of less than 5%, reliable data for 
Cd can be obtained, using the momentum-exchange method, for r > 5. That is, ten lattice spacings across 
the diameter of the cylinder are necessary to obtain reliable values of the force. This is consistent with the 
finding by Ladd [14]. In the range of 5 < r < 7, the stress-integration method produces larger fluctuations 
in the results than the momentum exchange method. For smaller radius, i.e., coarser resolutions, while both 
methods give poor results (due to insufficient resolution), the stress integration method yields much larger 
errors. 

Figure 7(b) compares Cd obtained from the methods of momentum exchange and the stress integration 
for Re = 10. The momentum-exchange method seems to gives a converged result at larger r (> 8). Based on 
the data for r > 8, an average value of Cd w 3.356 is obtained. In contrast, the stress-integration method has 
a larger fluctuation than the large r result from the momentum-exchange method even for r > 8. Averaging 
over the results for r > 8, the stress integration gives Cd & 3.319. The difference between converged results 
of two methods is about 1%. For r less than or around 5, the fluctuation in Cd from the stress-integration 
method is much larger than that in the momentum-exchange method. The conclusions from the comparisons 
in Fig. 7 are as follows: (i) both methods for force evaluation can give accurate results; (ii) the momentum- 
exchange method gives more consistent drag; and (iii) in the range of 10 < Re < 100, a resolution of ten 
lattice spacings across the diameter of the cylinder are needed in order to obtain consistent and reliable 
drag values. In other words, the lattice (grid) Reynolds number Re* (= U / v) should be less than 10 in the 
calculations. 

In the above results presented in Figs. 7(a) and 7(b), the center of the cylinder is placed on a lattice 
grid, thus the computational mesh is symmetric with respect to the geometry of the cylinder. To test the 
effect of the mesh symmetry on the accuracy of the force evaluation, the calculation of the flow' at Re = 10 
is repeated with different values of the cylinder center offset \ x in the x direction, or A„ in the y direction. 
The radius of the cylinder is deliberately chosen to be only 6.4 lattice grids. In order to preserve the mirror 
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Fig. 7. The drag coefficient for a uniform flow past a column of cylinders over a range of radius r. (a) Re = 100. The 
dashed line indicates the value of Cd — 1.24 obtained in Ref. [38]. (b) Re — 10. The dashed lines indicate the values of Co 
averaged over 4 largest radii. 


symmetry of the flow in the y direction, we use different boundary conditions for upper and lower boundaries 
(at y = ±77/2). For A„ = 0 while varying A x , we use the symmetric boundary conditions, which maintain 
the flow symmetry with respect to the center line in the x direction. For A* = 0 while varying A. y , we use 
the periodic boundary conditions at y = ±H/2, which are equivalent to the symmetric boundary conditions 
when A y = 0, but better reflect the flow symmetry when Ay ^ 0. The results of the drag coefficient Cd 
are presented in Table 3.2. The variation of Cd due to the change of the center of cylinder offset from a 
grid point is less than 1% when the cylinder diameter is only about 13 lattice spacings. The outcome is 
consistent with the expected truncation errors caused by mesh perturbation. We notice that the variation 
in Cd due to A x is about one order of magnitude smaller than that due to Ay. This is precisely because 
when Ay = 0 the mesh symmetry coincides with the flow symmetry in the y direction, and when A v y= 0 
the mesh symmetry is lost. This asymmetry due to Ay ^ 0 results in the change of the lift coefficient from 
0( 10~ 14 ) to O(10"“ 2 ), which is the same order of magnitude of the variation in Cd- It is our observation 
that the accuracy of the force evaluation schemes used here is dictated by that of the boundary conditions 
at the solid walls. The error due to symmetry of the computational mesh with respect to the geometry of 
an object is well bounded. This is also observed in other independent studies [23, 33]. 

Table 3.2 

The effect of symmetry of the computational mesh on the force evaluation for the steady uniform flow over a column of 
cylinders . The Reynolds number Re — 10 (r — 0.6,). the radius of the cylinder r — 6.4 (in the lattice unit of 6 X — }-), and 
H/r — 20. The variation of Co due to the change of the center of cylinder offset from a grid point is less than 1%. 


A* 

= 0, periodic boundary conditions at y 

= ±77/2 

A , 

0 0.2 

0.4 

0.6 

0.8 

Cd 

3.3661 3.3637 

3.3526 

3.3526 

3.3637 

Ay 

= 0, symmetric boundary conditions in y 

= ±77/2 

A, 

0 0.2 

0.4 

0.6 

0.8 

Cd 

3.3661 3.3666 

3.3646 

3.3667 

3.3692 


It is worth noting that the wall shear stress in the channel flow obtained by using the method of 
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momentum exchange has a relative error proportional to the resolution across the channel. For a resolution 
of 10 - 20 lattice spacings across the diameter considered here, the relative error in the drag appears, however, 
smaller than in the channel flow case. At Re = 100, with r > 10, the average value of the drag obtained by 
using the method of momentum exchange has a 1.7% relative error comparing with Fornberg’s data [38]. If 
the boundary layer thickness is estimated roughly to be 3 x 2r/\/Re as 6, there are only about six lattice 
spacings across the boundary layer over which the velocity profile changes substantially. Rased on the insight 
from the channel flow result, it is possible that the deviatoric shear stresses on the surface of the cylinder 
that are effectively incorporated in the method of momentum exchange suffer comparable levels of error as 
in the channel flow. The effective error cancellation over the entire surface of the body may have contributed 
to the good convergence behavior in the drag shown in Figs. 7(a) and 7(b). 

3.3. Flow over an asymmetrically placed circular cylinder in channel with vortex shedding. 

Schafer and Turek [39] reported a set of benchmark results for a laminar flow 7 over a circular cylinder of radius 
r that is asymmetrically placed inside a channel. In the present study, r = 12.8 is used and the center of 
the cylinder coincides with a grid point. The distance from the center of the cylinder to the upper wall 
arid lower wall is h+ — 4.2 r and h- = 4.0 r, respectively. This results in A + = 0.76 for the upper wall and 
A_ — 0.2 for the lower wall, respectively. The channel inlet has a parabolic profile and is placed at four 
radii upstream of the cylinder center according to the specification of the benchmark test [39] . This results 
in A = 0.2 for the inlet boundary. A zeroth-order extrapolation for f a is used at the exit boundary that is 
located 40 radii downstream of the cylinder center. Thus there are a total of 564 x 105 square lattices in the 
flow 7 field. For Re = 2rU/v = 100 based on the average inlet velocity U, the use of relaxation time r = 0.55 
requires 0 m 0.0651. 

At this Reynolds number, the flow becomes unsteady and periodic vortex shedding is observed. Fig- 
ures 8(a), 8(b), and 8(c), respectively, show time-dependent behaviors of the lift coefficient 


and the drag coefficient Cd [see Eq. (3.9)], and the pressure difference 


A P = 


P£ ~Pb 

Put'- ’ 


where pf arid pi, are the pressures at the front and the back of the cylinder, respectively, and po is the constant 
density imposed at the entrance. The data of Cl, Cd, and A P are compared with the benchmark results in 
Ref. [39]. We first note that the present numerical value of Strouhal number St = 2r/UT is about 0.3033, 
where T is the period of the lift curve. This agrees very well with the range of St values (0.2950 - 0.3050) 
given in Ref. [39]. We note that the difference in Ci,(t) between the momentum-exchange method and the 
surface stress-integration method is indiscernible graphically. For the drag coefficient CdU), it is interesting 
to note that although there is about 0.25% difference between the results given by the momentum-exchange 
method and the surface stress-integration method, both methods of force evaluation give two peaks in the 
Cr>(t) curves. Physically, these two peaks in the Cn(t) curve correspond to the existence of a weaker vortex 
and a stronger vortex alternately shed behind the cylinder. The difference in the strength of the vortices 
results from the difference: h+/r — 4.2 and h-/r = 4.0 in the passages between the cylinder and the channel 
walls. There is no report on the occurrence of these two peaks in Ref. [39]. Instead, a range of the maximum 
Cd (from 3.22 to 3.24) by different researchers was given. The present value of the higher peak is well within 
the range. It is interesting to note that both peaks of Co(t) obtained by the momentum-exchange method 
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Pig. 8. The 2D flow past a cylinder asymmetrically placed in a channel. The variations of the lift coefficient Cl, the 
drag coefficient Co, and the pressure difference AP as functions of time t (after an initial run time to) are compared with 
the benchmark results in Ref. [39]. At the lime to, the lift coefficient. Ci(t) attains its maximum value C“‘ ax . The dashed 
horizontal lines indicate the upper and lower bounds in Ref. [39]. The solid and dashed curves are the results obtained by rising 
momentum exchange and stress integration, respectively, (a) The lift coefficient Note that the results obtained by using 

the two methods are indistinguishable on the graph, (b) The drag coefficient Coit). (c) The pressure difference A P(t). The 
symbol x indicates the value of AP(to + T/2) given in Table 3.3, where T (fs 1296.5) is the period of Ciftt). 


are also within the range, as shown in Fig. 8(b). A further refined computation of the present problem using 
a multiblock procedure [40] with r = 40 in the fine grid region yield nearly the same results for C\o(t) and 
C L (t). 

We compile in Table 3.3 the values of Strouhal number St, maximum and minimum drag coefficient 
^- 7 max an( -] Cjj m , maximum and minimum lift coefficient and C™ m , and the pressure difference AP 
obtained by the LBE methods and other schemes of computational fluid dynamics given in Ref. [39]. The 
value of AP is measured at to +T/2, where to is the moment when (7 l(£) reaches its maximum value (7™ ax , 
and T is the periodicity of CT(£). For the LBE simulations, T is between 1296 and 1297 (in the lattice unit 
of St = 1). We use T = 1296.5 in the determination of the Strouhal number St. With a resolution much 
coarser than those used in Ref. [39], the LBE results are well within the bounds given in Ref. [39]. This 
clearly demonstrates the accuracy of the lattice Boltzmann method. 

3.4. Pressure-driven flow in a circular pipe. The steady-state flow field was obtained by using 
D3Q19 model with r = 0.52 [11]. Eq. (2.9) is used to evaluate the force on the boundary points along the 
circumference of the pipe over a distance of one lattice in the axial direction. The resulting axial force F x 
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Table 3.3 


Values of St, (7^ iax , Cjrf 111 , C^ iax , C£ mi , and 1 A P for the flow over a 2D cylinder asymmetrically placed in a channel. 
“ Momentum ” and “Stress” denote, respectively, the momentum- exchange method and the stress-integration method in the LBE 
calculations. The GFD results are the bounds in Ref. [39 j, which does not have data for Cn im and CV 11U . 


method 

St 

Qmax 


r»max 

rfmin 

A P 

Momentum 

0.3033 

3.2358 

3.1771 

1.0045 

-1.0347 

2.4914 

Stress 

0.3033 

3.2275 

3.1708 

1.0040 

-1.0340 

2.4914 

CFD 

0.2950 - 0.3050 

3.2200 - 3.2400 

— 

0.9900 - 1.0100 

— 

2.4600 - 2.5000 


is. equivalently, the force given by r w 2? rrS x , where r w is the wall shear stress and r is the pipe radius. For 
a fully developed flow inside a circular pipe, the exact fluid shear stress at the pipe wall is given by 

7~® xact (2 tit) - 7rr" . (3.10) 

ax 

We examine the normalized axial force, 


V = 


tv r- 


dp ■ 
dx 


(3.11) 


Figure 9 shows the normalized coefficient rj over a range of r: 3.5 - 23.5. Except for r < 5, rj is rather 
close to 1. It was noticed in Ref. [11] that the accuracy of LBE solution for the pipe flow is not as good 
as that for the two-dimensional channel flow due to the distribution of values of A around the pipe. The 
accuracy of the drag is dictated by the accuracy of the flow field if the force evaluation method is exact. 
For the pipe flow, the error in F x results from the inaccuracy in the flow field and the errors in the force 
evaluation scheme based on momentum exchange (as seen in the previous section for the two-dimensional 
channel flow case). For r > 5, the largest error in F x is about 3.5% and it occurs at r = 15.5. Again, 
there is no systematic error in F x . Given the complexity of the boundary in this three-dimensional flow, the 
results shown in Fig. 9 are satisfactory in the sense that it adds further credence to the momentum-exchange 
method for force evaluation. 



Fig. 9. The ratio r; between the tangential force F x on the pipe and its exact value (tit 2 V p) over a range of pipe radius r. 
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3.5. Steady uniform flow over a sphere. To limit the computational effort, a finite domain of 
—Hj 2 < y < Hj 2 and -H/2 < z < H/2, with H/r — 10 is used to compute the flow past a sphere of 
radius r (see Fig. 6). Two cases are considered: (a) the flow past a single sphere, and (b) the flow-' over a 
two-dimensional array of spheres (all located at x = 0) with the center of the spheres forming square lattices. 
In the former case, the boundary conditions at j v = 1 (y — H/2 corresponds to j v = 2) for f a 's are given by 
the following linear extrapolation 

faijx , 1 Jz) = 2 f a {jx,2,j z ) - fa(jx,3Jz) - (3.12) 


The velocity at j y — 2 is set as 


u(j x ,2,j z ) = u(j x ,S,j z ) . 


(3.13) 


Similar treatment is applied at y = H/2 and z = ±H/2. In the latter case, symmetry conditions are posed 
on f a ’s at j y = 1 by using the values of /„’ s at j v = 3 (see Ref. [28] for the two-dimensional case). At the 
inlet, a uniform velocity profile is imposed at j x = 1.5 (half way between the first and second lattices). The 
upstream boundary is located at 7.5 radii to the left of the sphere center in all simulations. 

For flow over a sphere, the drag coefficient, is often expressed as 


Cd 


24 


-4>, <i> 


F x 


(3.14) 
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where <j> accounts for the non-Stokesian effect of the drag. For two types of the boundary conditions at 
(: y — ±H/2 and z — ±H/2), 4> s denotes the non-Stokesian correction for the case where the symmetry 
conditions are imposed at ( y — ±H/2 and z = ±H/ 2) and 4<x> denotes the results for the case where the 
extrapolation for f a is used at ( y = ±H/2 and z = ±H/2) in order to simulate the unbounded flow. 

Figure 10(a) shows the non-Stokesian coefficient, </ :x , for r = 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 5.1, 5.2, 5.4, 5.6, 
and 5.8, for H/r = 10 at Re = 10. The relaxation time is r = 0.7. With this range of r, the number of 
the boundary nodes on the surface of the sphere increases roughly by a factor of (S.8/3) 2 « 3.74; the actual 
counts of the boundary nodes xr, gives a ratio 2370/546 = 4.35. The largest difference is 1.9% between r = 3.0 
and r = 3.2 that have the least resolution in the cases investigated. For a uniform flow over an unbounded 
sphere, an independent computation using a finite difference method based on the vorticity-stream function 
formulation with high resolution gives a drag coefficient <j> as 1.7986 at Re — 10. The largest difference 
between this result and the LBE results is 1.36% at r — 3.2. If the LBE data for the drag is averaged over 
the range of r, one obtains <p as 1.8086, which differs from 1.7986 by 0.54%. Hence, the LBE solutions with 


3.0 < r < 5.8 yield very consistent values for the drag force. Figure 10(b) shows the non-Stokesian correction 
factor <pg for ft uniform flow over a planar array of spheres for 3.0 < r < 5.8 and H/r = 10, at Re = 10. It 
is important to note that with the improvement of the surface resolution by a factor of 4.35, there is little 
systematic variation in < p s (r). The largest deviation from the average value, <p s as 1.963, is 1.1% at r = 5.0. It 
is clear that the LBE solution gives reliable fluid dynamic forces on a sphere at r as 3.5 for a moderate value 
of Re. The set of data for <j> s is inherently more consistent than that for </> x since the symmetry boundary- 
condition can be exactly specified at y = ±H/2 and 2 = ±H/2, while the extrapolation conditions given by 
Eqs. (3.12) and (3.13) do not guarantee the free stream condition at y = ±H/2 and 2 = ±H/2. Yet, both 


< {>00 and (p s exhibit remarkable self-consistency from coarse to not-so-coarse resolutions. 


4. Conclusions. Two methods for evaluating the fluid force in conjunction with the method of lattice 
Boltzmann equation for solving fluid flows involving curved geometry have been examined. The momentum- 
exchange method is very simple to implement. It is shown in the channel flow simulation that momentum- 
exchange method is not an exact method. The error in the wall shear stress is inversely proportional to 
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Fig. 10. Flcm past sphere. Variation of the non-Stokesian correction factor <p — — F x /§i:rU pv as a function of sphere 
radius r at Re — 10. The dashed lines are values of <f>(r) averaged over r. (a) The flow past a single sphere in an 'unbounded 
field (H/r — oo). (b) The flow past a planar array of spheres (H/r — 1 0 ) . 


the resolution. In two- and three-dimensional flows over a bluff body, it can give accurate drag values 
when there are at least ten lattice spacings across the body at Re ~ 100. The method of integrating the 
stresses on the surface of the body gives similar results when there is sufficient resolution but it exhibits 
much larger fluctuations than that in the method of momentum exchange when the resolution is limited. In 
addition, the stress-integration method requires considerably more efforts in implementing the extrapolation 
and integration on the body surface in comparison with the method of momentum exchange. 

It is interesting to note that the momentum-exchange method is perhaps superior to the stress-integration 
method because the former method is directly based on the distribution functions while the latter is derived 
from further processing of the distribution functions. In addition, the momentum-exchange method uses 
interpolations while the stress-integration method uses extrapolations. Often extrapolations are more noisy 
and unstable than interpolations. Even with a coarse resolution that does not yield very accurate local 
information, accurate force evaluation can be accomplished with the lattice Boltzmann method. Among the 
two force evaluation methods, the method of momentum exchange is recommended for force evaluation on 
curved boundaries for its simplicity, accuracy, and robustness. 
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Appendix A. LBE Models in Two and Three Dimensions. 

The nine- velocity (or 9-bit) LBE model on a two-dimensional square lattice, denoted as the D2Q9 model, 
has been widely used for simulations of two-dimensional flows. For three-dimensional flows, there are several 
cubic lattice models, such as the fifteen-velocity (D3Q15), nineteen-velocity (D3Q19), and twenty-seven- 
velocity (D3Q27) models, which have been used in the literature [10]. All these models have a rest particle 
(with zero velocity) in the discretized velocity set {e a \a = 0, 1, .... (b — 1)}. For athermal fluids, the 
equilibrium distributions for the D2Q9, D3Q15, D3Q19, and D3Q27 models are all of the following form 
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where w a is a weighting factor and e a is a discrete velocity, c = £ x /5 t is the unit speed, and S x and St are 
the lattice constant and the time step, respectively. The discrete velocities for the D2Q9 models are 


(0, 0), a = 0, 

(±1, 0) c, (0, ±1) c, a = 1, 3, 5, 7, 
(±1, ±1) c, a = 2, 4, 6, 8, 


and the values of the weighting factor w a are 


For the D3Q19 model, the discrete velocities are 

(0, 0), 

e a = \ (±1, 0, 0) c, (0, ±1, 0) c, (0, 0, ±1) c, 

(±1, ±1, 0) c, (0, ±1, ±1) c, (±1, 0, ±l)c, 


a = 0 , 
a — 1 — 6 , 
a — 7-18 , 


and the weighting factor w a is given by 


U'o 


i 
3 ' 

J_ 

18 ' 

1 

t 36 ’ 


a — 0, 
a — 16, 
a = 7-18. 


The discrete velocity sets {e a } for the D2Q9 and D3Q19 models are shown in Fig. 1. 
The density and velocity can be computed from , 

^=E/« = E/^ eQ) > 

a a 

p u = YL ea ^ a = Yl 

a a 

The speed of sound of the above LRE models is 


c, = 


VT 


and the equation of state is that of an ideal gas such that 


(A. 2) 



r j. 

a = 0, 



'U'a ~ S 

1 1 

1 9 ’ 

a = 1 , 1 

h 5, 7, 

(A. 3) 
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(A. 5) 


(A. 6a) 
(A. 6b) 


p = c;p. 


(A. 7) 


The viscosity of the fluid is 

" = <? s \ 

for the discrete velocity model of Eq, (1.2). It should be noted that the equilibrium distribution function 
is in fact a Taylor series expansion of the Maxwellian / ,0 > [6, 7]. This approximation of /// ql in algebraic 
form makes the LBE method valid only in the incompressible flow limit ujc — > 0. 
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Equation (1.2) is often discretized in space x and time t into the lattice Boltzmann equation 


fa( x i + e a S t , t + St) — fa(Xi,t) ----- --[ ) - f^\xi,t)] , (A. 8) 

r 

where r = A/5*. For this LBGK model [1, 2], the viscosity in the Navier-Stokes equation derived from the 
above lattice Boltzmann equation is 

v={t- c 2 g S t . (A. 9) 

The —1/2 correction in the above formula for u comes from the second-order derivatives of f a when f a (xi + 
e a 6t,t+6t) in Eq. (A. 8) is expanded in a Taylor series in u. This correction in v makes the lattice Boltzmann 
method formally a second-order method for solving incompressible flows [7]. Obviously, the physical and 
computational stabilities require that r > 1/2. 
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